clear all;

cd ..; cd ..; cd ..;
load('code_main_model\6_implied_w_r\computation_results\data.mat','V_backward','nmkt','pdf_B_forward',...
    'gl_backward','gk_backward','adjust_backward','r_guess','w_guess','psi','Z','chi','alpha','nu','GDP',...
    'nb','ntheta','nz','gridb','gridtheta','gridz','nb_pdf','ntheta_pdf','gridb_pdf','gridtheta_pdf');
cd 'code_figure_table\OA\Figure OA10';

r_end_index=min(find(r_guess(:,1)==0))-1;
w_end_index=min(find(w_guess(r_end_index,:,1,1,1)==0))-1;

n_years=11;

intensive(1:n_years,1)=0;
extensive(1:n_years,1)=0;
GE(1:n_years,1)=0;


temp_fin(1:n_years, 1:nmkt, 1:nz, 1:nb, 1:ntheta)=0.0;
l(1:n_years, 1:nmkt, 1:nz, 1:nb, 1:ntheta)=0.0;
temp_adjust(1:n_years, 1:nmkt, 1:nz, 1:nb, 1:ntheta)=0.0;

for i_t=1:n_years
    for i_mkt=1:nmkt
        for i_b=1:nb
            for i_theta=1:ntheta
                for i_z=1:nz
                    if (adjust_backward(i_t, i_mkt, i_z, i_b, i_theta)==1 || adjust_backward(i_t, i_mkt, i_z, i_b, i_theta)==2)
                        l(i_t, i_mkt, i_z, i_b, i_theta)=1.0;
                        temp_fin(i_t, i_mkt, i_z, i_b, i_theta)=0.0;
                    elseif (adjust_backward(i_t, i_mkt, i_z, i_b, i_theta)==3 || adjust_backward(i_t, i_mkt, i_z, i_b, i_theta)==4)
                        l(i_t, i_mkt, i_z, i_b, i_theta)=0.0;
                        temp_fin(i_t, i_mkt, i_z, i_b, i_theta)=0.0;    
                    elseif (adjust_backward(i_t, i_mkt, i_z, i_b, i_theta)==5 || adjust_backward(i_t, i_mkt, i_z, i_b, i_theta)==6)
                        l(i_t, i_mkt, i_z, i_b, i_theta)=0.0;
                        temp_fin(i_t, i_mkt, i_z, i_b, i_theta)=1.0;    
                    end
                end
            end
        end
    end

    
    for i_b=1:nb
        for i_theta=1:ntheta
            for i_z=1:nz
                if (adjust_backward(i_t, i_mkt, i_z, i_b, i_theta)==2 || adjust_backward(i_t, i_mkt, i_z, i_b, i_theta)==4 || adjust_backward(i_t, i_mkt, i_z, i_b, i_theta)==6)
                    temp_adjust(i_t, i_mkt, i_z, i_b, i_theta)=1.0;
                end
            end
        end
    end
end

for i_t=1:n_years
    
    r=r_guess(r_end_index,i_t);
    for i_mkt=1:nmkt
        w=w_guess(r_end_index,w_end_index,i_t,i_mkt);
        for i_b=1:nb_pdf
            for i_theta=1:ntheta_pdf
                for i_z=1:nz
                    if (pdf_B_forward(i_t,i_mkt,i_z,i_b,i_theta)~=0.0 || pdf_B_forward(1,i_mkt,i_z,i_b,i_theta)~=0.0)
    
                        k=lin_interpo2(nb, ntheta, gridb, gridtheta, squeeze(gk_backward(i_t,i_mkt,i_z,:,:)), gridb_pdf(i_b), gridtheta_pdf(i_theta));
                        labor=lin_interpo2(nb, ntheta, gridb, gridtheta, squeeze(gl_backward(i_t,i_mkt,i_z,:,:)), gridb_pdf(i_b), gridtheta_pdf(i_theta));
                        if ((k<1e-10) || (labor<1e-10))
                            k=0;
                            labor=0;
                        end
                        
                        k_prev=lin_interpo2(nb, ntheta, gridb, gridtheta, squeeze(gk_backward(1,i_mkt,i_z,:,:)), gridb_pdf(i_b), gridtheta_pdf(i_theta));
                        labor_prev=lin_interpo2(nb, ntheta, gridb, gridtheta, squeeze(gl_backward(1,i_mkt,i_z,:,:)), gridb_pdf(i_b), gridtheta_pdf(i_theta));
                        if ((k_prev<1e-10) || (labor_prev<1e-10))
                            k_prev=0;
                            labor_prev=0;
                        end
                        
                        y=0;
                        if ((k>0.0) && (labor>0.0))
                            y=Z(i_mkt)*gridz(i_z)*(k^alpha*labor^(1-alpha))^(1-nu);
                        end
                        
                        y_prev=0;
                        if ((k_prev>0.0) && (labor_prev>0.0))
                            y_prev=Z(i_mkt)*gridz(i_z)*(k_prev^alpha*labor_prev^(1-alpha))^(1-nu);
                        end
    
                        fin= lin_interpo2(nb, ntheta, gridb, gridtheta, squeeze(temp_fin(i_t,i_mkt,i_z,:,:)), gridb_pdf(i_b), gridtheta_pdf(i_theta))*psi(i_t,i_mkt);
                        fin_prev= lin_interpo2(nb, ntheta, gridb, gridtheta, squeeze(temp_fin(1,i_mkt,i_z,:,:)), gridb_pdf(i_b), gridtheta_pdf(i_theta))*psi(i_t,i_mkt);
                        
                        b=gridb_pdf(i_b);
                        theta=gridtheta_pdf(i_theta);
                        m=b*theta;
                        a=b*(1-theta);
                        if (k>m)
                            borrowing=k+psi(i_t,i_mkt)-m;
                        else
                            borrowing=0.0;
                        end
                        
                        if (k_prev>m)
                            borrowing_prev=k_prev+psi(1,i_mkt)-m;
                        else
                            borrowing_prev=0.0;
                        end
                             
                        production=(y-fin-borrowing*chi)*pdf_B_forward(i_t, i_mkt, i_z, i_b, i_theta);
                        production_prev=(y_prev-fin_prev-borrowing_prev*chi)*pdf_B_forward(1, i_mkt, i_z, i_b, i_theta);
                                           
                        if (borrowing_prev==0 && borrowing==0)
                            GE(i_t)=GE(i_t)+(production-production_prev);
                        elseif (borrowing_prev==0 && borrowing>0)
                            extensive(i_t)=extensive(i_t)+(production-production_prev);
                        elseif (borrowing_prev>0 && borrowing==0)
                            extensive(i_t)=extensive(i_t)+(production-production_prev);
                        elseif (borrowing_prev>0 && borrowing>0)
                            intensive(i_t)=intensive(i_t)+(production-production_prev);
                        end
    
                    end
                end
            end
        end
    end
end

GDP=sum(GDP,2);

GE=GE./GDP(1);
intensive=intensive./GDP(1);
extensive=extensive./GDP(1);

save('PB.mat','GE','intensive','extensive')